library(ggplot2)
library(reshape2)
library(Hmisc)
library(plyr)
library(latticeExtra)
library(MASS)
library(lme4)
library(lmerTest)
library(simpleboot)
library(gridExtra)
library(cowplot)
library(xtable)
library(ggridges)



toBePub <- read.csv("toBePub.csv")
relevantData = subset(toBePub, toBePub$enactmentYearPercentile>=0 & toBePub$enactmentYearPercentile<=100)
minYear <- aggregate(relevantData$Year, by= list(relevantData$Common_ID), FUN=min)
maxYear <- aggregate(relevantData$Year, by= list(relevantData$Common_ID), FUN=max)
meanTimeLasting <- mean(minYear$x - maxYear$x)
meanNruleChanges <- aggregate(relevantData$nRuleChanges, by= list(relevantData$Common_ID), FUN=sum)
meanNruleChanges <- mean(meanNruleChanges$x)




###### creating figure 1 #######

dat <- toBePub[toBePub$Country != "Spain",]
lev <- as.numeric(levels(as.factor(dat$Common_ID)))

# Overview of rule changes
dat$commons <- NA
for (cnt in 0:1){
  for (com in 1:10){
    dat$commons[dat$Common_ID==lev[10*cnt + com]] <- com
  }
}

tp <- dat[dat$enactmentYearPercentile >= 0 & dat$enactmentYearPercentile <= 100,]

new <- data.frame(commons=NA, country=NA, year=NA)
k <- 0
for(i in 1:nrow(tp)){
  if (tp$nRuleChanges[i] > 0){
    for (j in 1:tp$nRuleChanges[i]){
      k <- k + 1
      new[k,] <- c(tp$commons[i], tp$Country[i], tp$Year[i])
    }
  }
}

new$country <- factor(new$country, levels=c(1,3), labels=c("Netherlands", "United Kingdom"))
new$commons <- factor(new$commons)

g <- ggplot(new, aes(x=year, y=commons)) + theme_bw(base_size = 22)+ xlab("Year") + ylab("Commons") + theme(strip.background=element_rect(fill="#FFFFFF"))
plt <- g + geom_density_ridges(show.legend=F, scale = 0.9, bandwidth=7) + facet_grid(~ country, scales = "free_x")
plt


###### creating figure 2 #######

#this function creates a spline per common per category. The splines used in the figure are the mean of all
#the splines per common. Note, since there are different numbers of rule
#changes per common, this function has to be used to give every common equal weigth when drawing the figure.
mySmoother <- function (df, s=0.75)
{
  levelsCategories<-levels(df$Category)
  commons <- unique(as.factor(df$Common_ID))

  allCombinations <- expand.grid(categories = levelsCategories,Common_ID = commons)
  df$enactmentYearPercentile <- round(df$enactmentYearPercentile,0)

  smoothers <-
    lapply(1:length(allCombinations[,1]),
           function(n) ( loess(value~enactmentYearPercentile,
                               data = subset(df,df$Category == allCombinations[n,1] & df$Common_ID== allCombinations[n,2]
                                             & df$enactmentYearPercentile >= 0 & df$enactmentYearPercentile <= 100), span = s)
                        )
              )
  prediction <- data.frame()
  i <- 1
  for (j in levelsCategories)
  {
    indices <- seq(i,length(levelsCategories)*(length(commons)-1)+i,length(levelsCategories))
    smootherCategory <- smoothers[indices]
    newPrediction <-sapply(1:length(commons),function(n) (predict(smootherCategory[[n]],1:100)))
    newPrediction <- ifelse(newPrediction<0,0,newPrediction)
    nextCategory <- data.frame(enactmentYearPercentile = 1:100, smooth = rowMeans(newPrediction),Category=j)
    ifelse (i == 1, prediction <- nextCategory,prediction <- rbind(prediction,nextCategory))
    i <- i + 1
  }
  return(prediction)
}

bla <- melt(relevantData[,c(-3,-16,-17)], id.vars = c("Common_ID", "Year","enactmentYearPercentile","Country"), variable.name = c("Category"))

isOtherCategory <- bla$Category=="GeneralBureaucratic"|bla$Category=="GeneralResource"|bla$Category=="AppointmentBureaucratic"|bla$Category=="AppointmentResource"
bla$Category <- as.factor(ifelse(isOtherCategory,"other",as.character(bla$Category)))
bla$Category = factor(bla$Category,levels(bla$Category)[c(6,7,4,5,1,2,3)])
bla <- bla[order(bla$Category), ]

#####Plotting!!!


blaNL <- subset(bla,bla$Country=="Netherlands")
smoothedNL <- mySmoother(blaNL)
smoothedNL$country <- "Netherlands"
blaUK <- subset(bla,bla$Country=="United Kingdom")
smoothedUK <- mySmoother(blaUK)
smoothedUK$country <- "United Kingdom"
smoothed <- rbind(smoothedNL, smoothedUK)

labs <- c("Prohibition-Administrative", "Prohibition-Resources","Permission-Administrative", "Permission-Resources","Obligation-Administrative", "Obligation-Resources","Other")

g <- ggplot(smoothed, aes(enactmentYearPercentile,smooth,fill=Category)) +
  theme_bw(base_size = 22) + labs(x = "Standardised time",y="Average number of rule changes per time unit")
plt <- g + geom_area(position="stack")  + theme(legend.position = "top", legend.title = element_blank(), legend.spacing.x = unit(0.6,"mm")) + scale_fill_viridis_d(option="C", labels=labs) + facet_wrap(~country)

plt

#################################
### table 1
set.seed(123)
bootsize <- 1000
myBoot <- function(myModel)
{
  return(bootMer(myModel, fixef, nsim = bootsize, seed = 123,parallel="multicore",ncpus=4))
}
boot2stat <- function(myBoot)
{
  stat<-data.frame(cbind(beta=myBoot$t0,sigma=apply(myBoot$t,2,sd)))
  stat<-within(stat,{t<-beta/sigma;p<-sprintf("%.3f",2*pnorm(-abs(t)));t<-sprintf("%.3f",t);sigma<-sprintf("%.3f",sigma);beta<-sprintf("%.6f",beta)})[,c(1,2,4,3)]
  return(stat)
}


relevantData2 = subset(toBePub,toBePub$Country != "Spain")
relevantData2$Year = relevantData2$Year - min(relevantData2$Year) 

me1 <- lmer(nRuleChanges~ Year + I(Year*Year) + Country 
            + Year : Country + I(Year*Year) : Country + (1|Common_ID)
            ,data=relevantData2)
meBoot<-myBoot(me1) 
statMe <- boot2stat(meBoot)
rownames(statMe)<-c("$Intercept$","$Year$","$Year^2$", "UK","$Year$ X $UK$","$Year^2$ X $UK$")

print(xtable(statMe),sanitize.rownames.function=function(x) x,sanitize.colnames.function=function(x) x)


#################################
### table 2
me2 <- lmer(RatioBureaucratic~ Year + Country + nRuleChanges
            + Year : Country + (1|Common_ID)
            ,data=relevantData2)
meBoot<-myBoot(me2) 
statMe <- boot2stat(meBoot)
rownames(statMe)<-c("$Intercept$","$Year$","UK","N. of rule changes","$Year$ X $UK$")

print(xtable(statMe),sanitize.rownames.function=function(x) x,sanitize.colnames.function=function(x) x)


#################################
### table 3
me3 <- lmer(RatioSanction~ Year + Country + nRuleChanges
            + Year : Country + (1|Common_ID)
            ,data=relevantData2)
meBoot<-myBoot(me3) 
statMe <- boot2stat(meBoot)
rownames(statMe)<-c("$Intercept$","$Year$","UK","N. of rule changes","$Year$ X $UK$")

print(xtable(statMe),sanitize.rownames.function=function(x) x,sanitize.colnames.function=function(x) x)
